Survival Analysis on PBC Dataset
Kinga Bielicka, Kasper Kurzyński & Tomasz Warnke
22.01.2022
Wstęp
W niniejszej analizie wykorzystano zbiór danych pbc z pakietu Survival. Badany jest wpływ poszczególnych zmiennych, takich jak m.in. wiek, stadium choroby czy związane z badanym schorzeniem symptomy na przeżycie osób cierpiących na pierwotne stwardniające zapalenie dróg żółciowych. Pierwotne stwardniające zapalenie dróg żółciowych (ang. primary sclerosing cholangitis) jest chorobą autoimmunologiczną prowadzącą do wyniszczenia dróg żółciowych w wątrobie. Rozwój choroby jest powolny, lecz nieunikniony, prowadzący do marskości, a nawet wyniszczenia wątroby. Dane pochodzą z badania w Mayo Clinic przeprowadzonego w latach 1974-1984. Jednostką badania jest zbiorowość 424 pacjentów chorujących na PBC skierowanych na badanie pomiędzy 1974 a 1984 rokiem. W analizie przeżycia pod uwagę wzięto tylko pierwsze 312 przypadków, dla których dane są najbardziej kompletne. Zbiór pbc zawiera takie zmienne jak:
- Id – numer pacjenta
- Age – wiek w latach
- Albumin – serum albumin (g/dl), główne białka surowicy krwi, które są produkowane w wątrobie. Pełnią funkcję transportową dla hormonów, leków, kwasów tłuszczowych lub aminokwasów. Przybliżona norma dla osoby dorosłej to: 3,5 – 5,5 g/dl.
- Alk.phos – alkaline phosphatase (U/liter), fostaza alkaliczna to enzym występujący w całym organizmie. Jest to rodzaj białka w komórce, które działa jak katalizator i umożliwia zachodzenie pewnych procesów. Prawidłowy zakres fostazy alkaicznej w organizmie to 44 – 147 U/liter lub 30 – 120 U/liter według innych źródeł.
- Ascites – presence of ascites, wodobrzusze to nagromadzenie płynu w jamie brzusznej. Często występuje w wyniku marskości wątroby. Pacjenci posiadający wodobrzusze zostali oznaczeni jako „1”, natomiast pozostali jako „0”.
- Ast – aspartate aminotransferase (U/ml), test aminostransferazy asparaginianowej ma za zadanie dowiedzieć się czy wątroba jest uszkodzona. Prawidłowy wynik aktywności aminotransferazy asparaginianowej powinien mieścić się w granicach 5 – 40 IU/L. Przekroczona norma może wskazywać na ciężkie uszkodzenie wątroby.
- Bili – serum bilirunbin (mg/dl), bilirubina to barwnik żółciowy, który pochodzi z rozpadu krwinek czerwonych. Jej stężenie we krwi pozwala ocenić funkcjonowanie wątroby. Wzrost tego barwnika może być przyczyną żółtaczki. Bilirubina całkowita powinna mieścić się u osoby dorosłej w normach od 0,2 mg/dl do 1,1 mg/dl.
- Chol – serum cholesterol (mg/dl), Ilość cholesterolu w surowicy. Jego poziom u zdrowej osoby powinien wynosić poniżej 200 mg/dl.
- Copper – urine copper (ug/day), ilość miedzi w moczu, norma mieści się w przedziale 10 – 30 ug/day. Ilość miedzi w moczu bada się między innymi w związku z podejrzeniem marskości wątroby.
- Edema – obrzęk spowodowany nadmiarem płynu uwięzionego w tkankach organizmu. Często może być wynikiem marskości wątroby. Pacjenci niechorujący zostali oznaczeni jako „0”, pacjenci, u których schorzenie zostało wyleczone lub było nieleczone oznaczeni są jako „0,5”, natomiast pacjentów, u których edema występuje pomimo zastosowania terapii diuretykami oznaczono jako „1”.
- Hepato – obecność powiększonej wątroby (hepatomegalii), Pacjenci, u których zdiagnozowano powiększenie wątroby zostali oznaczeni jako „1”, natomiast pacjenci posiadający wątrobę o normalnych rozmiarach jako „0”.
- Platelet – ilość płytek krwi, inaczej trombocytów. Odgrywają kluczową rolę w procesie krzepnięcia. Za normę uważa się posiadanie 150 – 400 tys./mikrolitr.
- Protime – standardised blood clotting time, czas krzepnięcia krwi. Istnieją różne sposoby testowania czasu krzepnięcia krwi. Normą dla testu zastosowanego w badaniu Mayo jest przedział 9,5 – 11,8 sekundy.
- Sex – płeć, mężczyźni oznaczeni jako „m”, kobiety oznaczone jako „f”.
- Spiders – blood vessel malformations in the skin, malformacja naczyniowa to nienowotworowe zmiany naczyniowe, które spowodowane są nieprawidłowym rozwojem naczyń. Najczęściej pojawiają się na skórze i błonach śluzowych w obszarze głowy i szyi oraz w ośrodkowym układzie nerwowym. Pacjenci, u których została zdiagnozowana ta choroba zostali oznaczeni jako „1”, natomiast pozostali jako „0”.
- Stage – histologic stage of disease (needs biopsy) – Histologiczne stadium choroby I zapotrzebowanie na biopsję. Jego poziom został określony na skali od 1 do 4, gdzie 1 oznacza wczesne stadium, a 4 – późne.
- Status – status w punkcie końcowym. Został określony z pomocą trzech cyfr, gdzie „0” oznacza wartość cenzurowaną, „1” to pacjenci, którzy otrzymali przeszczep wątroby, natomiast „2” oznaczono osoby, które zmarły. Do celów projektowych pacjenci oznaczeni za pomocą „1” zostali zaklasyfikowani jako jednostki cenzurowane.
- Time – liczba dni pomiędzy rejestracją a śmiercią, transplantacją lub innym powodem przerwania badania.
- Trt – zmienna oznaczająca czy dany pacjent został poddany działaniu leku D-penicillmain lub placebo. „1” zostali oznaczeni pacjenci, którym został podany lek, natomiast „2” pacjenci, którzy otrzymali placebo.
- Trig – ilość trójglicerydów w organizmie. Trójglicerydy to tłuszcze proste i złożone, które są wykorzystywane jako budulec tkanki tłuszczowej i źródło energii. Ich podwyższony poziom może oznaczać między innymi stłuszczenie wątroby. Prawidłowy poziom trójglicerydów powinien wynosić poniżej 150 mg/dl.
Celem badania jest wyodrębnienie czynników wpływających na przeżywalność osób cierpiących na pierwotne stwardniające zapalenie dróg żółciowych.
Biblioteki wykorzystane do realizacji projektu:
- survival
- survminer
- survMisc
- dplyr
- gtools
- ggplot2
- corrplot
- gridExtra
- extrafont
- VIM
- mfp
- splines
- CoxR2
- ipred
- Rcpp
- RColorBrewer
- knitr
- kableExtra
- MASS
library(survival)
library(survminer)
library(survMisc)
library(dplyr)
library(gtools)
library(ggplot2)
library(corrplot)
library(gridExtra)
library(extrafont)
library(VIM)
library(extrafont)
library(mfp)
library(splines)
library(CoxR2)
library(ipred)
library(Rcpp)
library(RColorBrewer)
library(knitr)
library(kableExtra)
library(MASS)
data(pbc)
dane <- pbc
dane$cens <- ifelse(dane$status=="0" | dane$status=="1",0,1)
dane <- dane[1:312,]Prezentacja danych
Poniżej zaprezentowano 10 pierwszych i 10 ostatnich rekordów ze zbioru danych pbc. W zbiorze występują zmienne zarówno jakościowe, jak i ilościowe. Większość zmiennych jakościowych (wyjątek: sex) zakodowana została za pomocą cyfr “0” i “1” i na tym etapie pracy występują jako zmienne typu numerycznego. Do dalszej części projektu konieczne będzie ich przekształcenie.
head(dane,10)## id time status trt age sex ascites hepato spiders edema bili chol
## 1 1 400 2 1 58.76523 f 1 1 1 1.0 14.5 261
## 2 2 4500 0 1 56.44627 f 0 1 1 0.0 1.1 302
## 3 3 1012 2 1 70.07255 m 0 0 0 0.5 1.4 176
## 4 4 1925 2 1 54.74059 f 0 1 1 0.5 1.8 244
## 5 5 1504 1 2 38.10541 f 0 1 1 0.0 3.4 279
## 6 6 2503 2 2 66.25873 f 0 1 0 0.0 0.8 248
## 7 7 1832 0 2 55.53457 f 0 1 0 0.0 1.0 322
## 8 8 2466 2 2 53.05681 f 0 0 0 0.0 0.3 280
## 9 9 2400 2 1 42.50787 f 0 0 1 0.0 3.2 562
## 10 10 51 2 2 70.55989 f 1 0 1 1.0 12.6 200
## albumin copper alk.phos ast trig platelet protime stage cens
## 1 2.60 156 1718.0 137.95 172 190 12.2 4 1
## 2 4.14 54 7394.8 113.52 88 221 10.6 3 0
## 3 3.48 210 516.0 96.10 55 151 12.0 4 1
## 4 2.54 64 6121.8 60.63 92 183 10.3 4 1
## 5 3.53 143 671.0 113.15 72 136 10.9 3 0
## 6 3.98 50 944.0 93.00 63 NA 11.0 3 1
## 7 4.09 52 824.0 60.45 213 204 9.7 3 0
## 8 4.00 52 4651.2 28.38 189 373 11.0 3 1
## 9 3.08 79 2276.0 144.15 88 251 11.0 2 1
## 10 2.74 140 918.0 147.25 143 302 11.5 4 1tail(dane,10)## id time status trt age sex ascites hepato spiders edema bili chol
## 303 303 1250 0 2 60.65982 f 0 1 1 0 1.0 372
## 304 304 1230 0 1 35.53457 f 0 0 0 0 0.5 219
## 305 305 1216 0 2 43.06639 f 0 1 1 0 2.9 426
## 306 306 1216 0 2 56.39151 f 0 1 0 0 0.6 239
## 307 307 1149 0 2 30.57358 f 0 0 0 0 0.8 273
## 308 308 1153 0 1 61.18275 f 0 1 0 0 0.4 246
## 309 309 994 0 2 58.29979 f 0 0 0 0 0.4 260
## 310 310 939 0 1 62.33265 f 0 0 0 0 1.7 434
## 311 311 839 0 1 37.99863 f 0 0 0 0 2.0 247
## 312 312 788 0 2 33.15264 f 0 0 1 0 6.4 576
## albumin copper alk.phos ast trig platelet protime stage cens
## 303 3.25 108 1190 140 55 248 10.6 4 0
## 304 3.93 22 663 45 75 246 10.8 3 0
## 305 3.61 73 5184 288 144 275 10.6 3 0
## 306 3.45 31 1072 55 64 227 10.7 2 0
## 307 3.56 52 1282 130 59 344 10.5 2 0
## 308 3.58 24 797 91 113 288 10.4 2 0
## 309 2.75 41 1166 70 82 231 10.8 2 0
## 310 3.35 39 1713 171 100 234 10.2 2 0
## 311 3.16 69 1050 117 88 335 10.5 2 0
## 312 3.79 186 2115 136 149 200 10.8 2 0Ilość braków danych
Sprawdzone zostały także braki danych, których najwięcej jest w przypadku zmiennej dotyczącej poziomu cholesterolu (28 braków) oraz ilości trójglicerydów w organizmie (30 braków). Nieco mniejszą liczbą braków danych obarczone były takie zmienne jak ilość miedzi w moczu (2 braki) oraz ilość płytek krwi (4 braki).
NA_count <- colSums(is.na(dane))
NA_count## id time status trt age sex ascites hepato
## 0 0 0 0 0 0 0 0
## spiders edema bili chol albumin copper alk.phos ast
## 0 0 0 28 0 2 0 0
## trig platelet protime stage cens
## 30 4 0 0 0barplot(NA_count)
Imputacja danych i przekształcanie zmiennych
W celu pozbycia się braków danych wykorzystano funkcję kNN z pakietu VIM – imputację danych wykonano za pomocą metody k-najbliższych sąsiadów (parametr k ustawiono jako 5). Co więcej, zmienne wcześniej zakodowane jako zmienne ilościowe, przekształcono w zmienne jakościowe (status, trt, ascites, spiders, hepato, stage, edema, cens).
daneKNN <- kNN(dane, variable = c("chol", "copper", "trig", "platelet"), k = 5)
daneKNN <- daneKNN[,1:21,drop = F]
daneKNN$age <- round(daneKNN$age,0)
levels(daneKNN$sex) <- c(1,0)
daneKNN$status <- as.factor(daneKNN$status)
daneKNN$trt <- as.factor(daneKNN$trt)
daneKNN$ascites <- as.factor(daneKNN$ascites)
daneKNN$spiders <- as.factor(daneKNN$spiders)
daneKNN$hepato <- as.factor(daneKNN$hepato)
daneKNN$stage <- as.factor(daneKNN$stage)
daneKNN$edema <- as.factor(daneKNN$edema)
daneKNN$cens <- as.factor(daneKNN$cens)
str(daneKNN)## 'data.frame': 312 obs. of 21 variables:
## $ id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ time : int 400 4500 1012 1925 1504 2503 1832 2466 2400 51 ...
## $ status : Factor w/ 3 levels "0","1","2": 3 1 3 3 2 3 1 3 3 3 ...
## $ trt : Factor w/ 2 levels "1","2": 1 1 1 1 2 2 2 2 1 2 ...
## $ age : num 59 56 70 55 38 66 56 53 43 71 ...
## $ sex : Factor w/ 2 levels "1","0": 2 2 1 2 2 2 2 2 2 2 ...
## $ ascites : Factor w/ 2 levels "0","1": 2 1 1 1 1 1 1 1 1 2 ...
## $ hepato : Factor w/ 2 levels "0","1": 2 2 1 2 2 2 2 1 1 1 ...
## $ spiders : Factor w/ 2 levels "0","1": 2 2 1 2 2 1 1 1 2 2 ...
## $ edema : Factor w/ 3 levels "0","0.5","1": 3 1 2 2 1 1 1 1 1 3 ...
## $ bili : num 14.5 1.1 1.4 1.8 3.4 0.8 1 0.3 3.2 12.6 ...
## $ chol : int 261 302 176 244 279 248 322 280 562 200 ...
## $ albumin : num 2.6 4.14 3.48 2.54 3.53 3.98 4.09 4 3.08 2.74 ...
## $ copper : int 156 54 210 64 143 50 52 52 79 140 ...
## $ alk.phos: num 1718 7395 516 6122 671 ...
## $ ast : num 137.9 113.5 96.1 60.6 113.2 ...
## $ trig : int 172 88 55 92 72 63 213 189 88 143 ...
## $ platelet: int 190 221 151 183 136 233 204 373 251 302 ...
## $ protime : num 12.2 10.6 12 10.3 10.9 11 9.7 11 11 11.5 ...
## $ stage : Factor w/ 4 levels "1","2","3","4": 4 3 4 4 3 3 3 3 2 4 ...
## $ cens : Factor w/ 2 levels "0","1": 2 1 2 2 1 2 1 2 2 2 ...summary(daneKNN)## id time status trt age sex
## Min. : 1.00 Min. : 41 0:168 1:158 Min. :26.00 1: 36
## 1st Qu.: 78.75 1st Qu.:1191 1: 19 2:154 1st Qu.:42.00 0:276
## Median :156.50 Median :1840 2:125 Median :50.00
## Mean :156.50 Mean :2006 Mean :50.03
## 3rd Qu.:234.25 3rd Qu.:2697 3rd Qu.:57.00
## Max. :312.00 Max. :4556 Max. :78.00
## ascites hepato spiders edema bili chol
## 0:288 0:152 0:222 0 :263 Min. : 0.300 Min. : 120.0
## 1: 24 1:160 1: 90 0.5: 29 1st Qu.: 0.800 1st Qu.: 252.0
## 1 : 20 Median : 1.350 Median : 309.0
## Mean : 3.256 Mean : 364.8
## 3rd Qu.: 3.425 3rd Qu.: 394.2
## Max. :28.000 Max. :1775.0
## albumin copper alk.phos ast
## Min. :1.96 Min. : 4.00 Min. : 289.0 Min. : 26.35
## 1st Qu.:3.31 1st Qu.: 41.75 1st Qu.: 871.5 1st Qu.: 80.60
## Median :3.55 Median : 73.00 Median : 1259.0 Median :114.70
## Mean :3.52 Mean : 97.73 Mean : 1982.7 Mean :122.56
## 3rd Qu.:3.80 3rd Qu.:123.25 3rd Qu.: 1980.0 3rd Qu.:151.90
## Max. :4.64 Max. :588.00 Max. :13862.4 Max. :457.25
## trig platelet protime stage cens
## Min. : 33.00 Min. : 62.0 Min. : 9.00 1: 16 0:187
## 1st Qu.: 84.75 1st Qu.:200.0 1st Qu.:10.00 2: 67 1:125
## Median :107.50 Median :256.0 Median :10.60 3:120
## Mean :122.34 Mean :261.9 Mean :10.73 4:109
## 3rd Qu.:146.00 3rd Qu.:322.5 3rd Qu.:11.10
## Max. :598.00 Max. :563.0 Max. :17.10Graficzna prezentacja danych oraz wybrane statystyki opisowe
Rozkłady zmiennych ilościowych
a <- ggplot(data = daneKNN, aes(x = age)) + geom_histogram(alpha = 0.7, binwidth = 5, fill="darkorange") + theme_minimal() +
ggtitle('Age') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "N")
b <- ggplot(data = daneKNN, aes(x = bili)) + geom_histogram(alpha = 0.7, binwidth = 2, fill="mediumpurple3") + theme_minimal() +
ggtitle('Bilirunbin') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "N")
c <- ggplot(data = daneKNN, aes(x = chol)) + geom_histogram(alpha = 0.7, binwidth = 50, fill="seagreen4") + theme_minimal() +
ggtitle('Cholesterol') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "N")
d <- ggplot(data = daneKNN, aes(x = albumin)) + geom_histogram(alpha = 0.7, binwidth = 0.2, fill="dodgerblue3") + theme_minimal() +
ggtitle('Albumin') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "N")
a1 <- ggplot(data = daneKNN, aes(x = copper)) + geom_histogram(alpha = 0.7, binwidth = 30, fill="indianred4") + theme_minimal() +
ggtitle('Copper') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "N")
b1 <- ggplot(data = daneKNN, aes(x = alk.phos)) + geom_histogram(alpha = 0.7, binwidth = 400, fill="royalblue3") + theme_minimal() +
ggtitle('Alkaline phos.') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=20, family="serif")) +
labs(y = "N")
c1 <- ggplot(data = daneKNN, aes(x = ast)) + geom_histogram(alpha = 0.7, binwidth = 20, fill="goldenrod2") + theme_minimal() +
ggtitle('Aspartate amino.') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=20, family="serif")) +
labs(y = "N")
d1 <- ggplot(data = daneKNN, aes(x = trig)) + geom_histogram(alpha = 0.7, binwidth = 20, fill="cadetblue3") + theme_minimal() +
ggtitle('Triglycerides') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=20, family="serif")) +
labs(y = "N")
a2 <- ggplot(data = daneKNN, aes(x = platelet)) + geom_histogram(alpha = 0.7, binwidth = 25, fill="orangered3") + theme_minimal() +
ggtitle('Platelet') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=20, family="serif")) +
labs(y = "N")
b2 <- ggplot(data = daneKNN, aes(x = protime)) + geom_histogram(alpha = 0.7, binwidth = 0.5, fill="darkolivegreen4") + theme_minimal() +
ggtitle('Protime') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=20, family="serif")) +
labs(y = "N")
c2 <- ggplot(data = daneKNN, aes(x = time, fill = cens)) + geom_histogram(alpha = 0.7, binwidth = 200) + facet_grid(~cens) + theme_minimal() +
ggtitle('Time for cens & event') +
scale_fill_manual(values = c("coral3","slateblue3")) +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=20, family="serif")) +
labs(y = "N")Rozkład zmiennych age, bili, chol i albumin
W przypadku zmiennej dotyczącej wieku, najniższą wartością jest 26, najwyższą zaś 78. Średni wiek badanych osób to 50 lat i wokół tej wartości skupia się większość obserwacji. Poziom bilirubiny w normie waha się od 0,2 mg/dl do 1,1 mg/dl i w przypadku połowy obserwacji jej wartości nie przekraczają 1,35 mg/dl. Natomiast obserwacja, u której występuje najwyższa wartość bilirubiny to aż 28 mg/dl. Podobnie jest w przypadku cholesterolu, którego norma powinna wynosić mniej niż 200 mg/dl. Pomimo że rozkład zmiennej charakteryzuje się silną asymetrią prawostronną, to u większości jednostek cholesterol przekracza zakładaną normę, bowiem wartość pierwszego kwartyla to już 252 mg/dl. Wartość maksymalna przekracza normę niemal dziewięciokrotnie (1775 mg/dl). Zmienna albumin posiada rozkład lewostronny, choć w jego przypadku asymetria nie jest duża.
grid.arrange(a,b,c,d)
Rozkład zmiennych copper, alk.phos., asp i trig
Zmienne widoczne na poniższych wykresach charakterysują sie asymetrą prawostronną. Najmniejsza asymetria prawostronna widoczna jest dla zmiennej aspartate amino..
grid.arrange(a1,b1,c1,d1)
Rozkład zmiennych platelet i protime
Z histogramu dotyczącego zmiennej platelet wynika że w przypadku większości pacjentów ilość płytek krwi (trombocytów) była w normie (150 – 400 tys./mikrolitr). Podobnie jest w przypadku zmiennej protime, której 50% wartości mieściło się w przedziale [10, 11,1], a normą czasu krzepnięcia krwi w teście zastosowanym w badaniu Mayo jest przedział [9,5, 11,8] sekundy.
grid.arrange(a2,b2)
Rozkład zmiennej time z podziałem na jednostki cenzurowane i takie, które doznały zdarzenia
Rozkłady zmiennych dotyczących zarówno cenzurowania jak i zdarzenia względem czasu trwania nie różnią się znacząco.
c2
Rozkłady zmiennych jakościowych
e <- daneKNN %>%
count(sex) %>%
ggplot() + geom_col(aes(x = sex, y = n, fill = sex), alpha = 0.7)+ geom_label(aes(x = sex, y = n, label = n)) + scale_fill_manual(values = c("navajowhite3","plum4")) +
theme_minimal() +
ggtitle('Number of patients by gender') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
f <- daneKNN %>%
count(trt) %>%
ggplot() + geom_col(aes(x = trt, y = n, fill = trt), alpha = 0.7)+ geom_label(aes(x = trt, y = n, label = n)) + scale_fill_manual(values = c("lightseagreen","cornsilk3")) +
theme_minimal() +
ggtitle('Number of patients by treatment') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
g <- daneKNN %>%
count(ascites) %>%
ggplot() + geom_col(aes(x = ascites, y = n, fill = ascites), alpha = 0.7)+ geom_label(aes(x = ascites, y = n, label = n)) + scale_fill_manual(values = c("darkseagreen4","coral3")) +
theme_minimal() +
ggtitle('Number of patients by ascites presence') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
h <- daneKNN %>%
count(hepato) %>%
ggplot() + geom_col(aes(x = hepato, y = n, fill = hepato), alpha = 0.7)+ geom_label(aes(x = hepato, y = n, label = n)) + scale_fill_manual(values = c("palegreen4", "indianred3")) +
theme_minimal() +
ggtitle('Number of patients by hepato presence') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
e1 <- daneKNN %>%
count(edema) %>%
ggplot() + geom_col(aes(x = edema, y = n, fill = edema), alpha = 0.7)+ geom_label(aes(x = edema, y = n, label = n)) + scale_fill_manual(values = c("springgreen3","tan3","tomato3")) +
theme_minimal() +
ggtitle('Number of patients by edema presence') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
f1 <- daneKNN %>%
count(stage) %>%
ggplot() + geom_col(aes(x = stage, y = n, fill = stage), alpha = 0.7)+ geom_label(aes(x = stage, y = n, label = n)) + scale_fill_manual(values = c("palegreen3", "skyblue3","salmon3","brown3")) +
theme_minimal() +
ggtitle('Number of patients by stage of disease') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
g1 <- daneKNN %>%
count(spiders) %>%
ggplot() + geom_col(aes(x = spiders, y = n, fill = spiders), alpha = 0.7)+ geom_label(aes(x = spiders, y = n, label = n)) + scale_fill_manual(values = c("forestgreen", "chocolate3")) +
theme_minimal() +
ggtitle('Number of patients spiders presence') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
h1 <- daneKNN %>%
count(cens) %>%
ggplot() + geom_col(aes(x = cens, y = n, fill = cens), alpha = 0.7)+ geom_label(aes(x = cens, y = n, label = n)) + scale_fill_manual(values = c("coral3","slateblue3")) +
theme_minimal() +
ggtitle('Number of patients by cens or event') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))Podział pacjentów ze względu na zmienną sex i trt
Jeśli chodzi o rozkłady zmiennych jakościowych, obserwuje się znacznie większą liczbę pacjentów płci żeńskiej (276 jednostek) w porównaniu do liczby pacjentów płci męskiej (zaledwie 36 jednostek). Niemal połowa badanych jednostek otrzymała placebo (154 osoby), 158 osób zostało poddanych działaniu leku D-penicillmain.
grid.arrange(e,f)
Podział pacjentów ze względu na zmienną ascites i hepato
Spośród 312 badanych jednostek tylko 24 osoby cierpiały na wodobrzusze, natomiast jeżeli chodzi o obecność powiększonej wątroby, choroba ta została wykryta już u 160 pacjentów.
grid.arrange(g,h)
Podział pacjentów ze względu na zmienną edema i stage
U większości badanych (263 jednostki) obrzęk spowodowany nadmiarem płynu uwięzionego w tkankach organizmu nie występuje. W przypadku 29 osób schorzenie zostało wyleczone lub było nieleczone, a u 20 osób edema występuje pomimo zastosowania terapii diuretykami. Najwięcej badanych znajduje się w późniejszych stadiach choroby: III (120 jednostek) oraz IV (109 jednostek). Najmniej jednostek znajduje się w I stadium choroby i jest to 16 osób spośród wszystkich badanych.
grid.arrange(e1,f1)
Podział pacjentów ze względu na zmienną spiders i cens
Stosunkowo więcej jednostek niż w przypadku wyżej wymienionych schorzeń oraz objawów współwystępujących charakteryzuje się występowaniem zmian naczyniowych, bowiem to już 90 badanych. 125 spośród badanych jednostek to obserwacje cenzurowane.
grid.arrange(g1,h1)
Rozkłady pokazujące wybrane zależności
z <- ggplot(daneKNN, aes(x=stage, y = age, fill=stage)) + geom_boxplot(size=1.2, alpha=0.5) +
scale_fill_manual(values = c("palegreen3", "skyblue3","salmon3","brown3")) + theme_minimal() +
ggtitle("Age by stage") +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "Age")
x <- ggplot(daneKNN, aes(x=stage, y = bili, fill=stage)) + geom_boxplot(size=1.2, alpha=0.5) +
scale_fill_manual(values = c("palegreen3", "skyblue3","salmon3","brown3")) + theme_minimal() +
ggtitle("Bilirunbin by stage") +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "Bilirunbin")
z1 <- ggplot(daneKNN, aes(x=stage, y = time, fill=stage)) + geom_boxplot(size=1.2, alpha=0.5) +
scale_fill_manual(values = c("palegreen3", "skyblue3","salmon3","brown3")) + theme_minimal() +
ggtitle("Time by stage") +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "Time") + facet_grid(~cens)
x1 <- ggplot(daneKNN, aes(x=ascites, y = bili, fill=ascites)) + geom_boxplot(size=1.2, alpha=0.5) +
scale_fill_manual(values = c("darkseagreen4","coral3")) + theme_minimal() +
ggtitle("Bilirunbin by ascites presence") +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "Bilirunbin")
x2 <- ggplot(daneKNN, aes(x=hepato, y = bili, fill=hepato)) + geom_boxplot(size=1.2, alpha=0.5) +
scale_fill_manual(values = c("palegreen4", "indianred3")) + theme_minimal() +
ggtitle("Bilirunbin by hepato presence") +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "Bilirunbin")
z2 <- ggplot(daneKNN, aes(x=hepato, y = copper, fill=hepato)) + geom_boxplot(size=1.2, alpha=0.5) +
scale_fill_manual(values = c("palegreen4", "indianred3")) + theme_minimal() +
ggtitle("Copper by hepato presence") +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "Copper")Wykresy pudełkowe przedstawiające zależności pomiędzy zmienną age i stage oraz bili i stage
Najmłodsze z badanych osób znajdują się raczej w I stadium choroby, najstarsze zaś w ostatnim IV. Najwyższymi wartościami poziomu barwnika żółciowego charakteryzują się osoby z ostatniego stadia choroby.
grid.arrange(z,x)
Wykresy pudełkowe przedstawiające zależności pomiędzy zmienną time i stage oraz bili i ascites
W przypadku liczby dni pomiędzy rejestracją a śmiercią, transplantacją lub innym powodem przerwania badania – im wyższe stadium choroby, tym krótszy czas obserwacji.
Osoby z wodobrzuszem charakteryzują się dużo wyższymi wartościami stężenia barwnika żółciowego we krwi.
grid.arrange(z1,x1)
Wykresy pudełkowe przedstawiające zależności pomiędzy zmienną copper i hepato oraz bili i hepato
Podobna zależność występuje w przypadku związku ilości miedzi w moczu, a występowaniem powiększenia wątroby. Analogicznie jest również w przypadku zależności pomiędzy stężeniem barwnika żółciowego, a występowaniem powiększenia wątroby.
grid.arrange(z2,x2)
Statystyki względem grup
Statystyki dla wieku względem płci
Jeśli chodzi o różnice w wieku pomiędzy grupami kobiet i mężczyzn, nie są one wysokie, w grupie mężczyzn znajduje się jednostka najstarsza, natomiast w grupie kobiet jednostka najmłodsza. Średnia różni się pomiędzy grupami o 7 lat.
tab1 <- daneKNN %>%
group_by(sex) %>%
summarise(N = n(), Mean = round(mean(age),2), Min = min(age), Max = max(age), Sd = round(sd(age),2))
tab1 %>%
kbl() %>%
kable_styling(full_width = T, html_font = "Cambria", fixed_thead = T, font_size = 16)| sex | N | Mean | Min | Max | Sd |
|---|---|---|---|---|---|
| 1 | 36 | 56.19 | 33 | 78 | 11.54 |
| 0 | 276 | 49.22 | 26 | 77 | 10.20 |
Statystyki dla bilirunbiny względem obecności edemy
Stężenie barwnika żółciowego jest średnio wyższe u osób, u których występuje obrzęk spowodowany nadmiarem płynu uwięzionego w tkankach organizmu. Co więcej, odchylenie standardowe wskazuje na to, że wartości stężenia barwnika żółciowego u osób, u których edema nie występuje, nie różnią się tak bardzo wewnątrz tej grupy, jak w przypadku grup, u których schorzenie zostało wyleczone lub było nieleczone czy też u których edema występuje pomimo zastosowania terapii diuretykami.Stężenie barwnika żółciowego jest średnio wyższe u osób, u których występuje obrzęk spowodowany nadmiarem płynu uwięzionego w tkankach organizmu. Co więcej, odchylenie standardowe wskazuje na to, że wartości stężenia barwnika żółciowego u osób, u których edema nie występuje, nie różnią się tak bardzo wewnątrz tej grupy, jak w przypadku grup, u których schorzenie zostało wyleczone lub było nieleczone czy też u których edema występuje pomimo zastosowania terapii diuretykami.
tab2 <- daneKNN %>%
group_by(edema) %>%
summarise(N = n(), Mean = round(mean(bili),2), Min = min(bili), Max = max(bili), Sd = round(sd(bili),2))
tab2 %>%
kbl() %>%
kable_styling(full_width = T, html_font = "Cambria", fixed_thead = T, font_size = 16)| edema | N | Mean | Min | Max | Sd |
|---|---|---|---|---|---|
| 0 | 263 | 2.52 | 0.3 | 25.5 | 3.35 |
| 0.5 | 29 | 5.76 | 0.6 | 28.0 | 7.24 |
| 1 | 20 | 9.26 | 0.8 | 22.5 | 7.00 |
Statystyki dla bilirunbiny względem stopnia rozwoju choroby
Średnie stężenie barwnika żółciowego jest największe u osób, które posiadają najwyższe stadium rozwoju choroby.
tab3 <- daneKNN %>%
group_by(stage) %>%
summarise(N = n(), Mean = round(mean(bili),2), Min = min(bili), Max = max(bili), Sd = round(sd(bili),2))
tab3 %>%
kbl() %>%
kable_styling(full_width = T, html_font = "Cambria", fixed_thead = T, font_size = 16)| stage | N | Mean | Min | Max | Sd |
|---|---|---|---|---|---|
| 1 | 16 | 1.07 | 0.5 | 6.0 | 1.33 |
| 2 | 67 | 2.13 | 0.3 | 25.5 | 3.82 |
| 3 | 120 | 2.90 | 0.3 | 28.0 | 4.34 |
| 4 | 109 | 4.66 | 0.6 | 24.5 | 5.05 |
Statystyki dla poziomu cholesterolu względem płci
Kobiety charakteryzują się stosunkowo nieco wyższym poziomem cholesterolu niż mężczyźni, choć może to wynikać z dysproporcji w liczebności obu grup, bowiem również wartości skrajne są dużo wyższe.
tab4 <- daneKNN %>%
group_by(sex) %>%
summarise(N = n(), Mean = round(mean(chol),2), Min = min(chol), Max = max(chol), Sd = round(sd(chol),2))
tab4 %>%
kbl() %>%
kable_styling(full_width = T, html_font = "Cambria", fixed_thead = T, font_size = 16)| sex | N | Mean | Min | Max | Sd |
|---|---|---|---|---|---|
| 1 | 36 | 357.61 | 151 | 1000 | 178.80 |
| 0 | 276 | 365.79 | 120 | 1775 | 228.08 |
Statystyki ilości trójglicerydów względem płci
Jeśli chodzi o ilość trójglicerydów względem płci, średnia jest wyższa w przypadku mężczyzn, mimo że np. maksymalna wartość dla mężczyzn wynosi 242 mg/dl, a u kobiet jest dwukrotnie wyższa – 598 mg/dl.
tab5 <- daneKNN %>%
group_by(sex) %>%
summarise(N = n(), Mean = round(mean(trig),2), Min = min(trig), Max = max(trig), Sd = round(sd(trig),2))
tab5 %>%
kbl() %>%
kable_styling(full_width = T, html_font = "Cambria", fixed_thead = T, font_size = 16)| sex | N | Mean | Min | Max | Sd |
|---|---|---|---|---|---|
| 1 | 36 | 132.56 | 49 | 242 | 51.68 |
| 0 | 276 | 121.00 | 33 | 598 | 64.12 |
Statystyki ilości trójglicerydów względem stopnia rozwoju choroby
Średnie stężenie trójglicerydów we krwi rośnie również wraz ze stadium rozwoju choroby, choć średnia w III stadium jest nieco wyższa niż w stadium IV.
tab6 <- daneKNN %>%
group_by(stage) %>%
summarise(N = n(), Mean = round(mean(trig),2), Min = min(trig), Max = max(trig), Sd = round(sd(trig),2))
tab6 %>%
kbl() %>%
kable_styling(full_width = T, html_font = "Cambria", fixed_thead = T, font_size = 16)| stage | N | Mean | Min | Max | Sd |
|---|---|---|---|---|---|
| 1 | 16 | 90.75 | 55 | 188 | 32.22 |
| 2 | 67 | 112.63 | 44 | 319 | 53.19 |
| 3 | 120 | 128.68 | 46 | 382 | 56.91 |
| 4 | 109 | 125.96 | 33 | 598 | 75.28 |
Wykres pokazujący korelacje pomiędzy zmiennymi ilościowymi
Najsilniejsze korelacje występują w przypadku takich par zmiennych jak: protime, copper (0,22); ast, copper (0,3); trig, copper (0,29); ast, chol (0,35); trig, chol (0,28) czy też protime, platelet (-0,22); protime, albumin (-0,23); copper, albumin (-0,27); ast, albumin (-0,22).
doKor <- daneKNN[,c(5, 12,13,14,15,16,17,18,19)]
corMat <- cor(doKor)
#corrplot(corMat,method="number", tl.col = 'black')
corrplot(corMat, method = 'color',addCoef.col = 'black', order = 'AOE', tl.col = 'black', col=brewer.pal(n=10, name="PuOr"))
Funkcja czasu trwania oraz skumulowanego hazardu na podstawie zbioru danych pbc.
Porównanie estymatorów Kaplana - Meiera i Nelsona - Aelena
Na wykresach zostały pokazane funkcje przeżycia dla wszystkich jednostek z wykorzystaniem estymatora Kaplana - Meiera oraz Nelsona - Aelena. Dodatkowo, do oszacowania funkcji przeżycia za pomocą estymatora KM zostały wykorzystane dwa rodzaje przedziałów ufności - liniowy i logarytmiczny. Na wykresach brak jest istotnych różnic pomiędzy porównywanymi estymatorami. Potwierdza to również średnia oraz mediana czasu trwania widoczna na podsumowaniu pod wykresami. W momencie w którym wykresy się urywają (po 4000 dni) pozostaje już tylko około 25% jednostek.
daneKNN$cens <- as.character(daneKNN$cens)
daneKNN$cens <- as.numeric(daneKNN$cens)
pbcKM <- survfit(Surv(time, cens) ~ 1, conf.type="plain", data=daneKNN) #KM
pbcKM1 <- survfit(Surv(time, cens) ~ 1, conf.type="log-log", data=daneKNN)
pbcNELS=survfit(Surv(time,cens) ~ 1, conf.type ="plain", type="fleming-harrington",data=daneKNN)
splots <- list()
splots[[1]] <- ggsurvplot(
fit = pbcKM,
xlab = "Days",
ylab = "Overall survival probability with Kaplan - Meier estimator and plain confidence interval")
splots[[2]] <- ggsurvplot(
fit = pbcKM1,
xlab = "Days",
ylab = "Overall survival probability with with Kaplan - Meier estimator and log confidence interval",
palette = "blue")
splots[[3]] <- ggsurvplot(
fit = pbcNELS,
xlab = "Days",
ylab = "Overall survival probability with Nelson - Aelen estimator",
palette = "limegreen")
arrange_ggsurvplots(splots, print = TRUE,
ncol = 3, nrow = 1)
summary(pbcKM)$table## records n.max n.start events rmean se(rmean) median 0.95LCL
## 312.0000 312.0000 312.0000 125.0000 2973.6112 101.8206 3395.0000 3086.0000
## 0.95UCL
## 3839.0000summary(pbcKM1)$table## records n.max n.start events rmean se(rmean) median 0.95LCL
## 312.0000 312.0000 312.0000 125.0000 2973.6112 101.8206 3395.0000 3086.0000
## 0.95UCL
## 3839.0000summary(pbcNELS)$table## records n.max n.start events rmean se(rmean) median 0.95LCL
## 312.0000 312.0000 312.0000 125.0000 2978.9327 102.2104 3395.0000 3086.0000
## 0.95UCL
## 3839.0000daneKNN$Age_cat=quantcut(daneKNN$age,q=4)
daneKNN$bili_cat=quantcut(daneKNN$bili,q=4)
daneKNN$chol_cat=quantcut(daneKNN$chol,q=4)
daneKNN$albumin_cat=quantcut(daneKNN$albumin,q=4)
daneKNN$copper_cat=quantcut(daneKNN$copper,q=4)
daneKNN$alk_cat=quantcut(daneKNN$alk.phos,q=4)
daneKNN$ast_cat=quantcut(daneKNN$ast,q=4)
daneKNN$trig_cat=quantcut(daneKNN$trig,q=4)
daneKNN$plat_cat=quantcut(daneKNN$platelet,q=4)
daneKNN$pro_cat=quantcut(daneKNN$protime,q=4)Rozkłady czasu trwania dla grup
Grupy w zmiennych jakościowych zostały wydorębnione na podstawie wariantów zakodowanych w zmienych, natomiast zmienne ilościowe zostały podzielone za pomocą funkcji quantcut na 4 przedziały kwartylowe i każdy przedział stanowił odrębny wariant nowopowstałej zmiennej.
Wykresy funkcji trwania dla zmiennej sex i trt
Dla zmiennej sex większe prawdopodobieństwo przeżycia niemal w każdym momencie trwania jest w grupie kobiet. Mediana czasu przeżycia dla mężczyzn wynosi 2386 dni, natomiast dla kobiet 3428 dni. Test Log-Rank wskazał na istotne różnice pomiędzy tymi grupami jednostek, więc rozkłady czasu trwania są różne (p-value = 0.039). Z kolei dla pacjentów podzielonych względem zmiennej trt, test Log-Rank nie wskazał istotnych różnic (p-value = 0.75).
groupSplots1 <- list()
groupSplots1[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ sex, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("navajowhite3","plum4"))
groupSplots1[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ trt, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("lightseagreen","cornsilk3"))
arrange_ggsurvplots(groupSplots1, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmsex <- survfit(Surv(time, cens) ~ sex, data = daneKNN)
summary(kmsex)$table## records n.max n.start events rmean se(rmean) median 0.95LCL 0.95UCL
## sex=1 36 36 36 22 2479.267 281.2595 2386 1297 NA
## sex=0 276 276 276 103 3048.694 109.5150 3428 3170 NAkmtrt <- survfit(Surv(time, cens) ~ trt, data = daneKNN)
summary(kmtrt)$table## records n.max n.start events rmean se(rmean) median 0.95LCL 0.95UCL
## trt=1 158 158 158 65 2949.313 141.7174 3282 2583 NA
## trt=2 154 154 154 60 3002.749 145.7727 3428 3090 NAWykresy funkcji trwania dla zmiennej ascites i spiders
Pacjenci, u których stwierdzono wodobrzusze (ascites) posiadają dużo niższe prawdopodobieństwo przeżycia w każdym momencie trwania. Istnieją również istotne różnice w grupach wydzielonych na podstawie tej zmiennej (p-value < 0.0001). Podobna sytuacja klaruje się w przypadku zmiennej spiders (p-value < 0.0001). Mediana czasu trwania dla pacjentów, u których wykryto malformacje naczyniową wynosi zaledwie 1847 dni.
groupSplots2 <- list()
groupSplots2[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ ascites, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("darkseagreen4","coral3"))
groupSplots2[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ spiders, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("forestgreen", "chocolate3"))
arrange_ggsurvplots(groupSplots2, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmascites <- survfit(Surv(time, cens) ~ ascites, data = daneKNN)
summary(kmascites)$table## records n.max n.start events rmean se(rmean) median 0.95LCL
## ascites=0 288 288 288 102 3165.2653 102.3598 3584 3282
## ascites=1 24 24 24 23 856.1944 201.0002 368 223
## 0.95UCL
## ascites=0 NA
## ascites=1 1191kmspiders <- survfit(Surv(time, cens) ~ spiders, data = daneKNN)
summary(kmspiders)$table## records n.max n.start events rmean se(rmean) median 0.95LCL
## spiders=0 222 222 222 73 3287.546 113.0711 3839 3282
## spiders=1 90 90 90 52 2153.799 188.2438 1847 1434
## 0.95UCL
## spiders=0 NA
## spiders=1 3244Wykresy funkcji trwania dla zmiennej hepato i edema
Pacjenci, u których wykryto powiększoną wątrobę (hepato) i obrzęk spowodowany nadmiarem płynu uwięzionego w tkankach organizmu (edema) posiadają dużo niższe prawdopodobieństwo przeżycia w każdym momencie trwania. Również średnia dni trwania w poszczególnych grupach jest dużo niższa dla osób zmagających się z tymi chorobami. Test Log-Rank wykazał na istotne różnice pomiędzy grupami, a ich rozkładami czasu trwania (p-value < 0.0001).
groupSplots3 <- list()
groupSplots3[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ hepato, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("palegreen4", "indianred3"))
groupSplots3[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ edema, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("springgreen3","tan3","tomato3"))
arrange_ggsurvplots(groupSplots3, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmhepato <- survfit(Surv(time, cens) ~ hepato, data = daneKNN)
summary(kmhepato)$table## records n.max n.start events rmean se(rmean) median 0.95LCL 0.95UCL
## hepato=0 152 152 152 37 3606.189 129.0008 NA 3584 NA
## hepato=1 160 160 160 88 2372.598 138.0432 2256 1741 3244kmedema <- survfit(Surv(time, cens) ~ edema, data = daneKNN)
summary(kmedema)$table## records n.max n.start events rmean se(rmean) median 0.95LCL
## edema=0 263 263 263 89 3231.814 104.2501 3584 3244
## edema=0.5 29 29 29 17 2269.608 350.2181 1576 1012
## edema=1 20 20 20 19 673.550 207.7994 299 179
## 0.95UCL
## edema=0 NA
## edema=0.5 NA
## edema=1 971Wykresy funkcji trwania dla zmiennej stage
Wraz ze wzrostem stadium choroby spada prawdopodobieństwo przeżycia jednostek. Średnia czasu trwania dla pacjentów zmniejsza się, im stadium rozwoju choroby jest większe. Test Log-Rank wykazał istotne różnice rozkładów czasu trwania pomiędzy badanymi grupami (p-value < 0.0001).
ggsurvplot(survfit(Surv(time, cens) ~ stage, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("palegreen3", "skyblue3","salmon3","brown3"))
kmstage <- survfit(Surv(time, cens) ~ stage, data = daneKNN)
summary(kmstage)$table## records n.max n.start events rmean se(rmean) median 0.95LCL 0.95UCL
## stage=1 16 16 16 1 4358.727 188.0922 NA NA NA
## stage=2 67 67 67 16 3636.000 179.4216 4079 3839 NA
## stage=3 120 120 120 43 3156.503 153.5608 3428 2847 NA
## stage=4 109 109 109 65 2101.587 173.8259 1682 1165 2796Wykresy funkcji trwania dla skategoryzowanej zmiennej trig i age
W przypadku zmiennej trig prawdopodobieństwo przetrwania jest większe w grupach, dla których ilość trójglicerydów w organizmie jest niższa. Analogicznie jest w przypadku wieku - im niższy przedział wiekowy, tym większe jest prawdopodobieństwo przeżycia w danej grupie. Test Log-Rank wykazał istotne różnice rozkładów czasu trwania pomiędzy badanymi grupami (p-value = 0.0024 dla zmiennej trig oraz p-value < 0.0001 dla zmiennej age).
groupSplots4 <- list()
groupSplots4[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ trig_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"))
groupSplots4[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ Age_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_gray(),
palette = c("wheat3","seagreen3", "slateblue2", "sienna3"))
arrange_ggsurvplots(groupSplots4, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmtrig <- survfit(Surv(time, cens) ~ trig_cat, data = daneKNN)
summary(kmtrig)$table## records n.max n.start events rmean se(rmean) median
## trig_cat=[33,84.8] 78 78 78 22 3515.323 181.6379 NA
## trig_cat=(84.8,108] 78 78 78 32 2968.504 195.0771 3282
## trig_cat=(108,146] 81 81 81 30 2941.081 222.9211 3170
## trig_cat=(146,598] 75 75 75 41 2468.393 203.3881 2466
## 0.95LCL 0.95UCL
## trig_cat=[33,84.8] 3428 NA
## trig_cat=(84.8,108] 2796 NA
## trig_cat=(108,146] 3090 NA
## trig_cat=(146,598] 1682 4079kmage <- survfit(Surv(time, cens) ~ Age_cat, data = daneKNN)
summary(kmage)$table## records n.max n.start events rmean se(rmean) median 0.95LCL
## Age_cat=[26,42] 79 79 79 19 3492.100 199.5027 NA 3428
## Age_cat=(42,50] 85 85 85 28 3335.602 184.0169 4191 3358
## Age_cat=(50,57] 76 76 76 36 2831.014 192.2936 3282 2466
## Age_cat=(57,78] 72 72 72 42 2190.896 206.2281 2090 1356
## 0.95UCL
## Age_cat=[26,42] NA
## Age_cat=(42,50] NA
## Age_cat=(50,57] NA
## Age_cat=(57,78] 3222Wykresy funkcji trwania dla skategoryzowanej zmiennej bili i chol
W przypadku zmiennej bili prawdopodobieństwo przetrwania jest większe w grupach, dla których stężenie barwnika żółciowego w organizmie jest niższe. Podobnie jest w przypadku cholesterolu - im niższa jest wartość zmiennej chol, tym większe jest prawdopodobieństwo przeżycia w danej grupie, z wyjątkiem grupy, dla której cholesterol jest najniższy - w tym wypadku funkcja trwania pokrywa się z pozostałymi. Test Log-Rank wykazał istotne różnice rozkładów czasu trwania pomiędzy badanymi grupami (p-value < 0.0001 dla zmiennej bili oraz p-value = 0.00067 dla zmiennej chol).
groupSplots5 <- list()
groupSplots5[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ bili_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"))
groupSplots5[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ chol_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_gray(),
palette = c("wheat3","seagreen3", "slateblue2", "sienna3"))
arrange_ggsurvplots(groupSplots5, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmbili <- survfit(Surv(time, cens) ~ bili_cat, data = daneKNN)
summary(kmbili)$table## records n.max n.start events rmean se(rmean) median
## bili_cat=[0.3,0.8] 90 90 90 14 3964.305 139.3599 NA
## bili_cat=(0.8,1.35] 66 66 66 14 3722.350 188.9378 NA
## bili_cat=(1.35,3.42] 78 78 78 39 2683.480 180.4503 3222
## bili_cat=(3.42,28] 78 78 78 58 1447.146 144.4648 1170
## 0.95LCL 0.95UCL
## bili_cat=[0.3,0.8] NA NA
## bili_cat=(0.8,1.35] 3090 NA
## bili_cat=(1.35,3.42] 2105 3445
## bili_cat=(3.42,28] 943 1427kmchol <- survfit(Surv(time, cens) ~ chol_cat, data = daneKNN)
summary(kmchol)$table## records n.max n.start events rmean se(rmean) median
## chol_cat=[120,252] 79 79 79 28 3074.415 220.9585 NA
## chol_cat=(252,309] 78 78 78 24 3465.298 170.6644 3762
## chol_cat=(309,394] 77 77 77 28 2947.261 220.5058 2847
## chol_cat=(394,1.78e+03] 78 78 78 45 2411.435 180.8935 2105
## 0.95LCL 0.95UCL
## chol_cat=[120,252] 3170 NA
## chol_cat=(252,309] 3395 NA
## chol_cat=(309,394] 2419 NA
## chol_cat=(394,1.78e+03] 1492 3574Wykresy funkcji trwania dla skategoryzowanej zmiennej albumin i copper
W przypadku skategoryzowanej zmiennej albumin, najniższe prawdopodobieństwo przeżycia jest w grupie osób posiadających niski poziom wartości tej zmiennej. Test Log-Rank wykazał istotne różnice pomiędzy grupami w zależności od ich rozkładów czasu trwania (p-value < 0.0001). Podobna sytuacja występuje w przypadku zmiennej copper, jeśli chodzi o różne rozkłady czasu trwania (p-value < 0.0001). Jednakże dla tej zmiennej, niskie wartości oznaczają większe prawdopodobieństwo przeżycia.
groupSplots6 <- list()
groupSplots6[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ albumin_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"))
groupSplots6[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ copper_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_gray(),
palette = c("wheat3","seagreen3", "slateblue2", "sienna3"))
arrange_ggsurvplots(groupSplots6, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmalbumin <- survfit(Surv(time, cens) ~ albumin_cat, data = daneKNN)
summary(kmalbumin)$table## records n.max n.start events rmean se(rmean) median
## albumin_cat=[1.96,3.31] 81 81 81 55 1543.119 131.0277 1217
## albumin_cat=(3.31,3.55] 76 76 76 30 2985.155 203.5664 3358
## albumin_cat=(3.55,3.8] 81 81 81 18 3691.679 171.3645 NA
## albumin_cat=(3.8,4.64] 74 74 74 22 3495.126 170.9366 3762
## 0.95LCL 0.95UCL
## albumin_cat=[1.96,3.31] 1000 2256
## albumin_cat=(3.31,3.55] 2583 NA
## albumin_cat=(3.55,3.8] 3853 NA
## albumin_cat=(3.8,4.64] 3574 NAkmcopper <- survfit(Surv(time, cens) ~ copper_cat, data = daneKNN)
summary(kmcopper)$table## records n.max n.start events rmean se(rmean) median
## copper_cat=[4,41.8] 78 78 78 11 4043.353 140.0845 NA
## copper_cat=(41.8,73] 81 81 81 25 3250.959 191.9466 3358
## copper_cat=(73,123] 75 75 75 32 2875.992 205.1588 3170
## copper_cat=(123,588] 78 78 78 57 1734.077 172.5699 1170
## 0.95LCL 0.95UCL
## copper_cat=[4,41.8] NA NA
## copper_cat=(41.8,73] 3090 NA
## copper_cat=(73,123] 2400 NA
## copper_cat=(123,588] 930 1827Wykresy funkcji trwania dla skategoryzowanej zmiennej alk i ast
Skategoryzowane zmienne alk i ast wykazały, że jednostki należące do poszczególnych grup wyodrębnionych z uwagi na wartości tych zmiennych posiadają różne rozkładu czasu trwania. P-value w teście Log-Rank wyniosło odpowiednio 0.002 i < 0.0001. Pacjenci, którzy posiadają mniejsze wartości zmiennej alk i ast posiadają większe prawdopodobieństwo przeżycia, niż pozostali.
groupSplots7 <- list()
groupSplots7[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ alk_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"))
groupSplots7[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ ast_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_gray(),
palette = c("wheat3","seagreen3", "slateblue2", "sienna3"))
arrange_ggsurvplots(groupSplots7, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmalk <- survfit(Surv(time, cens) ~ alk_cat, data = daneKNN)
summary(kmalk)$table## records n.max n.start events rmean se(rmean)
## alk_cat=[289,872] 78 78 78 18 3589.857 198.9357
## alk_cat=(872,1.26e+03] 78 78 78 28 2942.289 201.0198
## alk_cat=(1.26e+03,1.98e+03] 78 78 78 31 2954.372 206.8846
## alk_cat=(1.98e+03,1.39e+04] 78 78 78 48 2454.019 186.7535
## median 0.95LCL 0.95UCL
## alk_cat=[289,872] NA NA NA
## alk_cat=(872,1.26e+03] 3222 2540 NA
## alk_cat=(1.26e+03,1.98e+03] 3574 2769 NA
## alk_cat=(1.98e+03,1.39e+04] 2256 1444 3358kmast <- survfit(Surv(time, cens) ~ ast_cat, data = daneKNN)
summary(kmast)$table## records n.max n.start events rmean se(rmean) median
## ast_cat=[26.4,80.6] 79 79 79 18 3721.073 166.1474 NA
## ast_cat=(80.6,115] 78 78 78 24 3299.247 196.2312 3853
## ast_cat=(115,152] 79 79 79 37 2600.008 200.7322 2796
## ast_cat=(152,457] 76 76 76 46 2225.745 197.0123 1690
## 0.95LCL 0.95UCL
## ast_cat=[26.4,80.6] 3762 NA
## ast_cat=(80.6,115] 3282 NA
## ast_cat=(115,152] 2288 3584
## ast_cat=(152,457] 1165 3445Wykresy funkcji trwania dla skategoryzowanej zmiennej plat i protime
W przypadku zmiennej plat prawdopodobieństwo przeżycia jest większe w grupach, w których pacjenci charakteryzowali się wyższą ilością trombocytów we krwi - funkcje przetrwania dla trzech przedziałów o najwyższych wartościach się pokrywają w pewnym stopniu, natomiast grupa, w której pacjenci mieli najmniejszą ilość płytek krwi, ma najmniejsze prawdopodobieństwo przeżycia. Test Log-Rank wykazał istotne różnice rozkładów czasu trwania pomiędzy badanymi grupami (p-value < 0.0001). Analogiczna sytuacja występuje w przypadku zmiennej dotyczącej krzepnięcia krwi - im dłuższy czas krzepnięcia krwi, tym prawdopodobieństwo przeżycia jest mniejsze. Test Log-Rank wykazał istotne różnice rozkładów czasu trwania pomiędzy badanymi grupami (p-value < 0.0001).
groupSplots8 <- list()
groupSplots8[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ plat_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"))
groupSplots8[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ pro_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_gray(),
palette = c("wheat3","seagreen3", "slateblue2", "sienna3"))
arrange_ggsurvplots(groupSplots8, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmplat <- survfit(Surv(time, cens) ~ plat_cat, data = daneKNN)
summary(kmplat)$table## records n.max n.start events rmean se(rmean) median
## plat_cat=[62,200] 80 80 80 48 2218.458 180.8581 2105
## plat_cat=(200,256] 77 77 77 27 3151.754 206.5502 3358
## plat_cat=(256,322] 77 77 77 22 3255.308 204.4453 3762
## plat_cat=(322,563] 78 78 78 28 3282.398 186.8769 NA
## 0.95LCL 0.95UCL
## plat_cat=[62,200] 1487 3428
## plat_cat=(200,256] 2503 NA
## plat_cat=(256,322] 3395 NA
## plat_cat=(322,563] 2847 NAkmpro <- survfit(Surv(time, cens) ~ pro_cat, data = daneKNN)
summary(kmpro)$table## records n.max n.start events rmean se(rmean) median
## pro_cat=[9,10] 89 89 89 16 3634.561 205.1407 4079
## pro_cat=(10,10.6] 85 85 85 24 3506.679 166.2759 3853
## pro_cat=(10.6,11.1] 61 61 61 31 2801.825 215.3994 3086
## pro_cat=(11.1,17.1] 77 77 77 54 1925.368 195.4251 1191
## 0.95LCL 0.95UCL
## pro_cat=[9,10] 4079 NA
## pro_cat=(10,10.6] 3222 NA
## pro_cat=(10.6,11.1] 2105 NA
## pro_cat=(11.1,17.1] 853 3090Graficzna prezentacjia funkcji skumulowanego hazardu dla wybranych zmiennych
Zmienne sex oraz stage
Hazard wzrasta szybciej w grupie mężczyzn niż w grupie kobiet. W przypadku zmiennej stage, największą wartość hazardu mają pacjenci z najwyższym stadium rozwoju choroby.
groupHplots <- list()
groupHplots[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ sex, data = daneKNN),
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_gray(),
palette = c("navajowhite3","plum4"),
fun = "cumhaz")
groupHplots[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ stage, data = daneKNN),
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_minimal(),
palette = c("palegreen3", "skyblue3","salmon3","brown3"),
fun = "cumhaz")
arrange_ggsurvplots(groupHplots, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
Zmienne hepato oraz skategoryzowana zmienna trig
Pacjenci, u których zostało stwierdzone powiększenie wątroby (hepato) mają większy hazard niż pacjenci bez tej dolegliwości. Wraz ze wzrostem czasu, różnica hazardów powiększa się. Dla zmiennej skategoryzowanej trig hazard jest najwyższy u pacjentów, którzy posiadają największe wartości trójglicerydów, lecz różnice nie są tak bardzo widoczne.
groupHplots <- list()
groupHplots[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ hepato, data = daneKNN),
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_gray(),
palette = c("palegreen4", "indianred3"),
fun = "cumhaz")
groupHplots[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ trig_cat, data = daneKNN),
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"),
fun = "cumhaz")
arrange_ggsurvplots(groupHplots, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
Skategoryzowane zmienne bili i copper
Hazard jest największy u osób posiadających dużą ilość miedzi w moczu oraz u osób, które posiadają największe stężenie barwnika żółciowego. Pacjenci należący do kategorii z mniejszymi wartościami obydwu zmiennych mają też odpowiednio niższe wartości hazardu.
groupHplots <- list()
groupHplots[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ copper_cat, data = daneKNN),
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_gray(),
palette = c("wheat3","seagreen3", "slateblue2", "sienna3"),
fun = "cumhaz")
groupHplots[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ bili_cat, data = daneKNN),
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"),
fun = "cumhaz")
arrange_ggsurvplots(groupHplots, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
Model PH Coxa
Modele regresji dla wybranych fukncji zmiennej losowej T, mają za zadanie pokazać relację między analizowaną funkcją, definiujacą rozkład zmiennej losowej T, a zbiorem zmiennych charakteryzujących jednostki należące do badanej populacji. W analizie przeżycia najczęściej modelowanymi funkcjami są funkcja hazardu, funkcja trwania, a czasami także dystrybuanta. W modelu zaproponowanych przez Coxa [1972], funkcja hazardu bazowego szacowana jest nieparametrycznie dla momentów tk, w których wystąpiło zdarzenie, czyli jest to model semiparametryczny. Model Coxa nazywany jest modelem proporcjonalnych hazardów (PH), co oznacza, że dla dwóch jednostek z wektorami zmiennych x i x*, iloraz ich hazardów nie zależy od czasu t.
- Zbudowanie modelu PH Coxa będzie zawierało następujące zagadnienia:
- interpretacja uzyskanych oszacowań parametrów modelu,
- weryfikacja założeń modelu (proporcjonalność hazardów),
- wybór postaci funkcji zmiennych objaśniających,
- selekcja zmiennych do modelu,
- reszty w modelu,
- wybór najlepszego modelu,
- ocena dopasowania modelu.
Estymacja modelu
Globalny test Walda (p=<2e-16) wskazuje, że przynajmniej 1 zmienna objaśniająca w modelu statystycznie różni sie od 0, czyli jest istotna statystycznie. Zmienne: sex, ascites, spiders, chol, alk.phos, ast, trig oraz platelet nie sa różne od 0, bo p-value > 0,05.
Cox <- coxph(Surv(time,cens)~age+sex+ascites+hepato+spiders+bili+chol+albumin+copper+alk.phos+ast+trig+platelet+protime, data = daneKNN)
summary(Cox) ## Call:
## coxph(formula = Surv(time, cens) ~ age + sex + ascites + hepato +
## spiders + bili + chol + albumin + copper + alk.phos + ast +
## trig + platelet + protime, data = daneKNN)
##
## n= 312, number of events= 125
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 2.706e-02 1.027e+00 1.033e-02 2.619 0.008815 **
## sex0 -1.019e-01 9.031e-01 2.972e-01 -0.343 0.731787
## ascites1 3.428e-01 1.409e+00 3.411e-01 1.005 0.314994
## hepato1 4.667e-01 1.595e+00 2.257e-01 2.068 0.038625 *
## spiders1 1.863e-01 1.205e+00 2.222e-01 0.838 0.401818
## bili 8.917e-02 1.093e+00 2.174e-02 4.102 4.1e-05 ***
## chol 1.604e-04 1.000e+00 3.945e-04 0.407 0.684291
## albumin -1.047e+00 3.511e-01 2.608e-01 -4.013 6.0e-05 ***
## copper 2.949e-03 1.003e+00 1.150e-03 2.564 0.010349 *
## alk.phos -2.074e-05 1.000e+00 3.810e-05 -0.544 0.586184
## ast 3.255e-03 1.003e+00 1.731e-03 1.880 0.060082 .
## trig -1.136e-03 9.989e-01 1.253e-03 -0.906 0.364814
## platelet -3.698e-04 9.996e-01 1.118e-03 -0.331 0.740890
## protime 3.053e-01 1.357e+00 8.253e-02 3.700 0.000216 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0274 0.9733 1.0068 1.0484
## sex0 0.9031 1.1072 0.5044 1.6172
## ascites1 1.4088 0.7098 0.7219 2.7494
## hepato1 1.5947 0.6271 1.0247 2.4818
## spiders1 1.2048 0.8300 0.7794 1.8623
## bili 1.0933 0.9147 1.0477 1.1409
## chol 1.0002 0.9998 0.9994 1.0009
## albumin 0.3511 2.8482 0.2106 0.5854
## copper 1.0030 0.9971 1.0007 1.0052
## alk.phos 1.0000 1.0000 0.9999 1.0001
## ast 1.0033 0.9968 0.9999 1.0067
## trig 0.9989 1.0011 0.9964 1.0013
## platelet 0.9996 1.0004 0.9974 1.0018
## protime 1.3571 0.7369 1.1544 1.5953
##
## Concordance= 0.851 (se = 0.019 )
## Likelihood ratio test= 187.6 on 14 df, p=<2e-16
## Wald test = 210.4 on 14 df, p=<2e-16
## Score (logrank) test = 317.7 on 14 df, p=<2e-16Zastosowanie funkcji krokowej, która poszukuje modelu najbardziej zoptymalizowanego ze wzgledu na kryterium akaike AIC (tu - metoda wsteczna)
W wyniku tej procedury otrzymano model z 7 zmiennymi objasniajacymi: age,hepato,bili,albumin,copper,ast,protime. Według kryterium AIC zmienna ast mimo p-value jest równe 0,06 powinna znalezc sie w modelu. Pozostałe zmienne na poziomie istotnosci 0,05 sa statystycznie rozne od 0 (czyli istotne), ponieważ p-value < 0,05.
stepAIC(Cox,direction="backward")## Start: AIC=1120.33
## Surv(time, cens) ~ age + sex + ascites + hepato + spiders + bili +
## chol + albumin + copper + alk.phos + ast + trig + platelet +
## protime
##
## Df AIC
## - platelet 1 1118.4
## - sex 1 1118.5
## - chol 1 1118.5
## - alk.phos 1 1118.6
## - spiders 1 1119.0
## - trig 1 1119.2
## - ascites 1 1119.3
## <none> 1120.3
## - ast 1 1121.7
## - hepato 1 1122.7
## - copper 1 1124.4
## - age 1 1125.2
## - protime 1 1129.6
## - bili 1 1132.2
## - albumin 1 1134.0
##
## Step: AIC=1118.44
## Surv(time, cens) ~ age + sex + ascites + hepato + spiders + bili +
## chol + albumin + copper + alk.phos + ast + trig + protime
##
## Df AIC
## - chol 1 1116.5
## - sex 1 1116.6
## - alk.phos 1 1116.8
## - spiders 1 1117.2
## - trig 1 1117.5
## - ascites 1 1117.5
## <none> 1118.4
## - ast 1 1120.6
## - hepato 1 1121.2
## - copper 1 1122.4
## - age 1 1123.5
## - protime 1 1127.9
## - bili 1 1130.4
## - albumin 1 1132.3
##
## Step: AIC=1116.53
## Surv(time, cens) ~ age + sex + ascites + hepato + spiders + bili +
## albumin + copper + alk.phos + ast + trig + protime
##
## Df AIC
## - sex 1 1114.7
## - alk.phos 1 1114.9
## - spiders 1 1115.3
## - trig 1 1115.5
## - ascites 1 1115.6
## <none> 1116.5
## - ast 1 1119.0
## - hepato 1 1119.3
## - copper 1 1120.4
## - age 1 1121.5
## - protime 1 1125.9
## - bili 1 1130.2
## - albumin 1 1130.3
##
## Step: AIC=1114.7
## Surv(time, cens) ~ age + ascites + hepato + spiders + bili +
## albumin + copper + alk.phos + ast + trig + protime
##
## Df AIC
## - alk.phos 1 1113.1
## - spiders 1 1113.4
## - trig 1 1113.6
## - ascites 1 1113.9
## <none> 1114.7
## - ast 1 1117.5
## - hepato 1 1117.6
## - copper 1 1120.2
## - age 1 1121.0
## - protime 1 1124.0
## - bili 1 1128.2
## - albumin 1 1128.4
##
## Step: AIC=1113.14
## Surv(time, cens) ~ age + ascites + hepato + spiders + bili +
## albumin + copper + ast + trig + protime
##
## Df AIC
## - trig 1 1112.1
## - spiders 1 1112.2
## - ascites 1 1112.7
## <none> 1113.1
## - hepato 1 1116.0
## - ast 1 1116.1
## - copper 1 1118.2
## - age 1 1120.6
## - protime 1 1122.1
## - albumin 1 1126.4
## - bili 1 1126.5
##
## Step: AIC=1112.11
## Surv(time, cens) ~ age + ascites + hepato + spiders + bili +
## albumin + copper + ast + protime
##
## Df AIC
## - ascites 1 1111.0
## - spiders 1 1111.2
## <none> 1112.1
## - hepato 1 1114.9
## - ast 1 1115.7
## - copper 1 1117.2
## - age 1 1120.9
## - protime 1 1121.9
## - bili 1 1125.2
## - albumin 1 1125.2
##
## Step: AIC=1111.05
## Surv(time, cens) ~ age + hepato + spiders + bili + albumin +
## copper + ast + protime
##
## Df AIC
## - spiders 1 1110.4
## <none> 1111.0
## - hepato 1 1113.6
## - ast 1 1114.2
## - copper 1 1118.5
## - protime 1 1121.7
## - age 1 1121.9
## - bili 1 1126.2
## - albumin 1 1128.5
##
## Step: AIC=1110.39
## Surv(time, cens) ~ age + hepato + bili + albumin + copper + ast +
## protime
##
## Df AIC
## <none> 1110.4
## - ast 1 1113.4
## - hepato 1 1114.0
## - copper 1 1119.1
## - age 1 1120.7
## - protime 1 1122.4
## - bili 1 1127.5
## - albumin 1 1129.1## Call:
## coxph(formula = Surv(time, cens) ~ age + hepato + bili + albumin +
## copper + ast + protime, data = daneKNN)
##
## coef exp(coef) se(coef) z p
## age 0.0326504 1.0331893 0.0092752 3.520 0.000431
## hepato1 0.5080490 1.6620454 0.2174273 2.337 0.019458
## bili 0.0854208 1.0891753 0.0177598 4.810 1.51e-06
## albumin -1.1054157 0.3310732 0.2348614 -4.707 2.52e-06
## copper 0.0032776 1.0032830 0.0009424 3.478 0.000505
## ast 0.0037099 1.0037168 0.0015580 2.381 0.017255
## protime 0.3220291 1.3799250 0.0778588 4.136 3.53e-05
##
## Likelihood ratio test=183.5 on 7 df, p=< 2.2e-16
## n= 312, number of events= 125Estymacja nowego modelu Coxa po odrzuceniu zmiennych nieistotnych statystycznie
- INTERPRETACJA PARAMETROW:
- jeżeli wiek wzrośnie o rok, hazard (ryzyko zgonu pod warunkiem, że jednostka dożyje do danego momentu) rośnie o 3,3%, ceteris paribus,
- pacjenci z obecnością hepatomegalii (powiększonej wątroby) mają hazard o 66,2% wyższy niż pacjenci z hepato = 0, ceteris paribus,
- wraz ze wzrostem bili o 1, hazard rośnie o 8,9%, ceteris paribus,
- wraz ze wzrostem albumin o 1, hazard maleje o 67%, ceteris paribus,
- wraz ze wzrostem copper o 1, hazard rośnie o 0,3%, ceteris paribus,
- wraz ze wzrostem ast o 1, hazard rośnie o 0,4%, ceteris paribus,
- wraz ze wzrostem protime o 1, hazard rośnie o 38,0%, ceteris paribus.
Cox1 <- coxph(Surv(time,cens)~age+hepato+bili+albumin+copper+ast+protime, data = daneKNN)
summary(Cox1) ## Call:
## coxph(formula = Surv(time, cens) ~ age + hepato + bili + albumin +
## copper + ast + protime, data = daneKNN)
##
## n= 312, number of events= 125
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.0326504 1.0331893 0.0092752 3.520 0.000431 ***
## hepato1 0.5080490 1.6620454 0.2174273 2.337 0.019458 *
## bili 0.0854208 1.0891753 0.0177598 4.810 1.51e-06 ***
## albumin -1.1054157 0.3310732 0.2348614 -4.707 2.52e-06 ***
## copper 0.0032776 1.0032830 0.0009424 3.478 0.000505 ***
## ast 0.0037099 1.0037168 0.0015580 2.381 0.017255 *
## protime 0.3220291 1.3799250 0.0778588 4.136 3.53e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0332 0.9679 1.0146 1.0521
## hepato1 1.6620 0.6017 1.0853 2.5452
## bili 1.0892 0.9181 1.0519 1.1278
## albumin 0.3311 3.0205 0.2089 0.5246
## copper 1.0033 0.9967 1.0014 1.0051
## ast 1.0037 0.9963 1.0007 1.0068
## protime 1.3799 0.7247 1.1846 1.6074
##
## Concordance= 0.845 (se = 0.019 )
## Likelihood ratio test= 183.5 on 7 df, p=<2e-16
## Wald test = 199.3 on 7 df, p=<2e-16
## Score (logrank) test = 276.9 on 7 df, p=<2e-16Weryfikacja założenia proporcjonalności
Kluczowym założeniem, które musi być spełnione, by model prawidłowo odzwierciedlał wpływ zmiennych objaśniających na rozkład czasu trwania, jest założenie proporcjonalności hazardów. Spełnienie tego założenia powinno być zweryfikowane dla każdej ze zmiennych objaśniających w modelu. Jeżeli p > 0,05 to zostaly spełnione zalożenia proporcjonalności. Dla zmiennych: bili oraz protime te założenia nie zostały spełnione, więc należy usunąć te zmienne z modelu.
Prop <- cox.zph(Cox1, transform = 'km')
print(Prop)## chisq df p
## age 0.3590 1 0.5491
## hepato 0.8706 1 0.3508
## bili 7.3198 1 0.0068
## albumin 1.1252 1 0.2888
## copper 0.0145 1 0.9043
## ast 1.9955 1 0.1578
## protime 5.2823 1 0.0215
## GLOBAL 16.9008 7 0.0180par(mfrow=c(4,2))
plot(Prop) 
Estymacja modelu z wykorzystaniem zmiennych istotnych statystycznie oraz spełniających założenie o proporcjonalności
- INTERPRETACJA PARAMETROW:
- jeżeli wiek wzrośnie o rok, hazard (ryzyko zgonu pod warunkiem, że jednostka dożyje do danego momentu) rośnie o 3,4%, ceteris paribus,
- pacjenci z obecnością hepatomegalii (powiększonej wątroby) mają hazard o 102% wyższy niz pacjenci z hepato = 0, ceteris paribus,
- wraz ze wzrostem albumin o 1, hazard maleje o 70%, ceteris paribus,
- wraz ze wzrostem copper o 1, hazard rośnie o 0,5%, ceteris paribus,
- wraz ze wzrostem ast o 1, hazard rośnie o 0,5%, ceteris paribus.
Cox2 <- coxph(Surv(time,cens)~age+hepato+albumin+copper+ast, data = daneKNN)
summary(Cox2) ## Call:
## coxph(formula = Surv(time, cens) ~ age + hepato + albumin + copper +
## ast, data = daneKNN)
##
## n= 312, number of events= 125
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.0330674 1.0336202 0.0088058 3.755 0.000173 ***
## hepato1 0.7046633 2.0231653 0.2112452 3.336 0.000851 ***
## albumin -1.2051939 0.2996339 0.2336910 -5.157 2.51e-07 ***
## copper 0.0046068 1.0046175 0.0008553 5.386 7.19e-08 ***
## ast 0.0047647 1.0047761 0.0012846 3.709 0.000208 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0336 0.9675 1.0159 1.0516
## hepato1 2.0232 0.4943 1.3373 3.0609
## albumin 0.2996 3.3374 0.1895 0.4737
## copper 1.0046 0.9954 1.0029 1.0063
## ast 1.0048 0.9952 1.0022 1.0073
##
## Concordance= 0.824 (se = 0.02 )
## Likelihood ratio test= 140.4 on 5 df, p=<2e-16
## Wald test = 162.9 on 5 df, p=<2e-16
## Score (logrank) test = 174.9 on 5 df, p=<2e-16Graficzna weryfikacja założenia proporcjonalności
Metoda skalowanych reszt Schoenfelda Ocenę tego, czy założenie proporcjonalności jest spełnione dokonuje się wzrokowo poprzez sprawdzenie czy punkty na wykresie układają się wzdłuż linii poziomej. Wyraźna zaobserwowana tendencja w przebiegu tej krzywej może sugerować brak proporcjonalności zmiennej.
Ocena na podstawie wykresy wydaje się być problematyczna, ale można zauważyć, że zmienne bill i protime nie spełniają założeń proporcjonalności.
reszty <- resid(Cox1, type="scaledsch")
Time <- as.numeric(rownames(reszty))
zmienne <- names(daneKNN[,c(5,8,11,13,14,16,19)])
par(mfrow = c(3,3))
for (i in 1:7) {
plot(log(Time), reszty[,i], xlab="ln(Czas)", main=zmienne[i],
ylab="Skalowane reszty Schoenfelda", pch=20, cex=0.7)
lines(smooth.spline(log(Time), reszty[,i] ), lwd=3 )
}
Identyfikacja jednostek odstających
Ważnym aspektem ewaluacji modeli Coxa jest badanie jednostek odstających, tj. tych które:
- mają nietypową konfigurację zmiennych objaśniających,
- wywierają niepożądany wpływ na oszacowania parametrów,
- mają niepożądany wpływ na dopasowanie modelu [Hosmer i inni 2008].
Są to punkty na wykresach oddalone od pozostałych. Jednostki o nietypowej konfiguracji zmiennych objaśniających. Mogą w nieproporcjonalny sposób wpływać na oszacowania parametrów strukturalnych modelu. Do identyifikacji wartości nietypowych wykorzystuje się reszty odchyleń. Za wartości nietypowe można uznać te jednostki, dla których reszty odchyleń przyjmują wartości +- 3 lub więcej.
deviance <- residuals(Cox2,type="deviance")
s <- Cox2$linear.predictors
plot(s,deviance,xlab="Liniowy predyktor",ylab="Reszty odchylen",cex=0.5, pch=20)
abline(h=c(3,-3),lty=3)
daneKNN$deviance <- deviance
c1 <- which(daneKNN$deviance < c(-3))Jednostki wpływowe
W przypadku wysokich wartości przyrostu oceny parametru należy rozważyć możliwość usunięcia danej jednostki z próby. Do usunięcia z próby zakwalifikowały się jednostki o numerze 48, 166, 231, 233 i 253. Są one odbiegające i wpływowe.
dfb <- residuals(Cox2,type="dfbeta")
n <- dim(dfb)[1]
obs.nr <- c(1:n)
par(mfrow = c(2,3))
for (j in 1:5) {
plot(obs.nr,dfb[,j],xlab="Numer jednostki",ylab="Przyrost oceny parametru",
main=zmienne[j])
}
a1 <- which(abs(dfb[,1])>(0.004)) #usunac te jednostki
a2 <- which(abs(dfb[,2])>(0.06))
a3 <- which(abs(dfb[,3])>(0.1))
a4 <- which(abs(dfb[,4])>(0.0004))
a5 <- which(abs(dfb[,5])>(0.001))
c <- sort(unique(c(a1,a2,a3,a4,a5,c1)))
daneKNN_1 <- daneKNN[-c,] #zredukowany zbior
Oszacowanie modelu na zredukowanej bazie danych
- INTERPRETACJA PARAMETROW:
- jeżeli wiek wzrośnie o 1 rok, hazard (ryzyko zgonu pod warunkiem, że jednostka dożyje do danego momentu) rośnie o 4,9%, ceteris paribus,
- pacjenci z obecnością hepatomegalii (powiększonej wątroby) mają hazard o 100,5% wyższy niż pacjenci z hepato = 0, ceteris paribus,
- wraz ze wzrostem albumin o 1, hazard maleje o 76%, ceteris paribus,
- wraz ze wzrostem copper o 1, hazard rośnie o 0,6%, ceteris paribus,
- wraz ze wzrostem ast o 1, hazard rośnie o 0,8%, ceteris paribus.
Cox3 <- coxph(Surv(time,cens)~age+hepato+albumin+copper+ast, data = daneKNN_1)
summary(Cox3)## Call:
## coxph(formula = Surv(time, cens) ~ age + hepato + albumin + copper +
## ast, data = daneKNN_1)
##
## n= 307, number of events= 124
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.0486020 1.0498024 0.0098801 4.919 8.69e-07 ***
## hepato1 0.6956171 2.0049458 0.2086548 3.334 0.000857 ***
## albumin -1.4274764 0.2399136 0.2545430 -5.608 2.05e-08 ***
## copper 0.0064491 1.0064699 0.0009382 6.874 6.25e-12 ***
## ast 0.0082022 1.0082359 0.0016792 4.885 1.04e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0498 0.9526 1.0297 1.0703
## hepato1 2.0049 0.4988 1.3320 3.0179
## albumin 0.2399 4.1682 0.1457 0.3951
## copper 1.0065 0.9936 1.0046 1.0083
## ast 1.0082 0.9918 1.0049 1.0116
##
## Concordance= 0.836 (se = 0.019 )
## Likelihood ratio test= 180.3 on 5 df, p=<2e-16
## Wald test = 187.2 on 5 df, p=<2e-16
## Score (logrank) test = 217.2 on 5 df, p=<2e-16METODA RESZT MARTYNGALOWYCH - badamy czy wystepuje liniowosc
Jeżeli zmienna objaśniająca 𝑋 jest ilościowa, to wprowadzenie jej do modelu Coxa jest równoważne przyjęciu założenia, że łączy ją związek liniowy z logarytmem hazardu. Jedna z najcześciej spotykanych metod graficznych jest metoda reszt martyngałowych. Polega na wyznaczeniu reszt martyngałowych dla modelu bez zmiennych objaśniających i zestawieniu ich z badaną zmienną objaśniającą na wykresie korelacyjnym.
Reszty martyngałowe dla zmiennej age
Jeżeli funkcja ta ma kształt zbliżony do liniowego to potwierdza to liniowość relacji między zmienną objaśniającą a logarytmem hazardu, w przeciwnym wypadku kształt funkcji może być wskazówką co do wyboru odpowiedniej metody transformacji zmiennej.
age - Potwierdzone założenie o liniowości
zmn1 <- daneKNN_1$age
lab1 <- "Age (years)"
reszty1 <- resid(coxph(Surv(time,cens)~1,data=daneKNN_1),type="martingale")
plot(zmn1, reszty1, xlab=lab1,ylab="Martingale residuals",cex=0.6)
lines(lowess(zmn1, reszty1,delta=1),lwd=2)
Reszty martyngałowe dla zmiennej albumin
albumin - Potwierdzone założenie o liniowości
zmn2 <- daneKNN_1$albumin
lab2 <- "albumin"
reszty2 <- resid(coxph(Surv(time,cens)~1,data=daneKNN_1),type="martingale")
plot(zmn2, reszty2, xlab=lab2,ylab="Reszty martynga?owe",cex=0.6)
lines(lowess(zmn2, reszty2,delta=1),lwd=2)
Reszty martyngałowe dla zmiennej copper
copper - Brak potwierdzenia założenia o liniowości
zmn3 <- daneKNN_1$copper
lab3 <- "copper"
reszty3 <- resid(coxph(Surv(time,cens)~1,data=daneKNN_1),type="martingale")
plot(zmn3, reszty3, xlab=lab3,ylab="Reszty martynga?owe",cex=0.6)
lines(lowess(zmn3, reszty3,delta=1),lwd=2)
Reszty martyngałowe dla zmiennej ast
ast - Potwierdzone założenie o liniowości
zmn4 <- daneKNN_1$ast
lab4 <- "ast"
reszty4 <- resid(coxph(Surv(time,cens)~1,data=daneKNN_1),type="martingale")
plot(zmn4, reszty4, xlab=lab4,ylab="Reszty martynga?owe",cex=0.6)
lines(lowess(zmn4, reszty4,delta=1),lwd=2)
Transformacja zmiennej copper
W przypadku podejrzenia braku liniowości związku, stosowane są różne metody transformacji zmiennej objaśniającej. Najpopularniejsze z nich, to: dychotomizacja zmiennej objaśniającej (wyznaczenie optymalnego punktu odcięcia)
Dychotomizacja zmiennej ciągłej to zastąpienie jej przez zmienną zero-jedynkową taką, że zmienna ta przyjmuje wartość 0, gdy zmienna ciągła przyjmuje wartości równe lub mniejsze niż punkt odcięcia, oraz wartość 1, gdy zmienna ciągła przyjmuje wartości większe niż punkt odcięcia.
Wartość kryterium informacyjnego Akaike po dychotomizacji okazała się dużo gorsza niż przed tym zabiegiem.
Cox4 <- coxph(Surv(time, cens)~zmn3,data=daneKNN_1)
A1 <- round(AIC(Cox3),2) #akaike z modelu coxa
A1## [1] 1093.72cutpoint <- cutp(Cox4)$zmn3 ## optymalny punkt odciecia
c <- (zmn3>=cutpoint$zmn3[1])*1+0 #jezeli zmienna jest >= punkotwi odciecia to przypisuje 1 w przeciwnym wypadku 0
Cox5 <- coxph(Surv(time, cens)~c,data=daneKNN_1)
A2 <- round(AIC(Cox5),2)
A2## [1] 1201.05Zastosowanie wielomianów ułamkowych
Metoda polega na poszukiwaniu dla zmiennej 𝑋, funkcji 𝑔(𝑥), która najlepiej odzwierciedla prawdziwą relację zmiennej 𝑋 i logarytmu funkcji hazardu.
Wartość kryterium informacyjnego Akaike ponownie jest wyższa w porównaniu do zmiennej copper przed transformacją.
m3 <- mfp(Surv(time,cens)~ fp(zmn3, df = 4, select = 0.05),family=cox, data=daneKNN_1) ## wielomiany ulamkowe
A3 <- round(AIC(m3))
A3## [1] 1182Metoda wielomianowej funkcji sklejanej (splajn)
int <- quantile(na.omit(zmn3),probs=c(0.05,0.275,.5,.725,.95))
Cox6 <- coxph(Surv(time,cens)~ns(zmn3,knots=int),data=daneKNN_1)
pred <- predict(Cox6,type="terms",se.fit=TRUE, terms=1)
plot(na.omit(zmn3),exp(pred$fit),type="n",xlab=lab3,ylim=c(0,2.5),
ylab="Hazard względny")
lines(smooth.spline(na.omit(zmn3),exp(pred$fit+1.96*pred$se.fit)),lty=2)
lines(smooth.spline(na.omit(zmn3),exp(pred$fit-1.96*pred$se.fit)),lty=2)
lines(smooth.spline(na.omit(zmn3),exp(pred$fit)),lty=1)
abline(h=1,lty=3)
legend('topright',2,c("splajn","95% przedzial ufnosci"), lty=1:2, box.lty=0)
A4 <- round(AIC(Cox6),2)
A4## [1] 1187.47Podsumowanie transformacji zmiennej copper
- A1 - model bez tranformacji
- A2 - po dychotomizacji
- A3 - wielomiany ułamkowe
- A4 - splajn
Mimo wszystko najlepszym modelem jest model bez transformacji zmiennej copper. Interpretacja parametrów w tym modelu przedstawia się w następujący sposób:
- jezeli wiek wzrosnie o 1 rok, hazard (ryzyko zgonu pod warunkiem, że jednostka dożyje do danego momentu) rośnie o 5,0%, ceteris paribus,
- pacjenci z obecnością hepatomegalii (powiększonej wątroby) mają hazard o 100,5% wyższy niz pacjenci z hepato = 0, ceteris paribus,
- wraz ze wzrostem albumin o 1, hazard maleje o 76%, ceteris paribus,
- wraz ze wzrostem copper o 1, hazard rośnie o 0,6%, ceteris paribus,
- wraz ze wzrostem ast o 1, hazard rośnie o 0,8%, ceteris paribus.
rbind(A1,A2,A3,A4)## [,1]
## A1 1093.72
## A2 1201.05
## A3 1182.00
## A4 1187.47Cox3 <- coxph(Surv(time,cens)~age+hepato+albumin+copper+ast, data = daneKNN_1)
summary(Cox3)## Call:
## coxph(formula = Surv(time, cens) ~ age + hepato + albumin + copper +
## ast, data = daneKNN_1)
##
## n= 307, number of events= 124
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.0486020 1.0498024 0.0098801 4.919 8.69e-07 ***
## hepato1 0.6956171 2.0049458 0.2086548 3.334 0.000857 ***
## albumin -1.4274764 0.2399136 0.2545430 -5.608 2.05e-08 ***
## copper 0.0064491 1.0064699 0.0009382 6.874 6.25e-12 ***
## ast 0.0082022 1.0082359 0.0016792 4.885 1.04e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0498 0.9526 1.0297 1.0703
## hepato1 2.0049 0.4988 1.3320 3.0179
## albumin 0.2399 4.1682 0.1457 0.3951
## copper 1.0065 0.9936 1.0046 1.0083
## ast 1.0082 0.9918 1.0049 1.0116
##
## Concordance= 0.836 (se = 0.019 )
## Likelihood ratio test= 180.3 on 5 df, p=<2e-16
## Wald test = 187.2 on 5 df, p=<2e-16
## Score (logrank) test = 217.2 on 5 df, p=<2e-16Miary dopasowania
Poniżej ukazano miary dopasowania R2, (obliczone na podstawie statystyki cząstkowej ilorazu wiarygodności w modelu Coxa) Najlepiej dopasowanym modelem (o największej wartości tej miary = 0,77) jest model Cox3 otrzymany metodą krokową, po odrzuceniu zmiennych niespełniajacych założeń proporcjonalnosci oraz jednostek odstających i wpływowych. Co ciekawe, taką samą wartość miary R2 (0,77) obliczono dla modelu Cox1 otrzymanego metodą krokową przed sprawdzeniem zalożeń proporcjonalnosci oraz przed odrzuceniem jednostek odstających i wpływowych
r2_1 <- coxr2(Cox1) #model otrzymany metoda krokowa przed sprawdzeniem zalozen proporcjonalnosci
r2_2 <- coxr2(Cox2) #model po odrzuceniu zmiennych niespelniajacych zalozen proporcjonalnosci
r2_3 <- coxr2(Cox3) #model po odrzuceniu jednostek odstajacych i wplywowych (uznany za najlepszy)
round(rbind(r2_1$rsq, r2_2$rsq, r2_3$rsq),2)## rsq
## [1,] 0.77
## [2,] 0.67
## [3,] 0.77